Matrix-valued Boltzmann equation for the non-integrable Hubbard chain 
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The standard Fermi-Hubbard chain becomes non-integrable by adding to the nearest neighbor 
hopping additional longer range hopping amplitudes. We assume that the quartic interaction is 
weak and investigate numerically the dynamics of the chain on the level of the Boltzmann type 
kinetic equation. Only the spatially homogeneous case is considered. We observe that the huge 
degeneracy of stationary states in case of nearest neighbor hopping is lost and the convergence to 
the thermal Fermi-Dirac distribution is restored. The convergence to equilibrium is exponentially 
fast. However for small n.n.n. hopping amplitudes one has a rapid relaxation towards the manifold 
of quasi-stationary states and slow relaxation to the final equilibrium state. 



I. INTRODUCTION 

The most widely known quantum chains are integrable, 
in the sense that they have a large number of local conser- 
vation laws. Eigenfunctions can be determined through 
the Bethe ansatz and there are special relations for scat- 
tering amplitudes, to mention only a few characteristics, 
see [1, 2] for further details. Obviously, dynamical prop- 
erties depend sensitively on the integrable structure. For 
example, such chains have a large Drude weight generi- 
cally, signaling ballistic transport but still leaving room 
for a diffusive component [3]. There has been consider- 
able efforts to understand what happens as one moves 
away from integrability [4, 5]. In our contribution we 
study the case where integrability is lost by adding cou- 
plings beyond the nearest neighbor ones. But we will 
stay in the regime where kinetic theory remains applica- 
ble. More than by other methods, we arrive at detailed 
information on how non-integrability becomes manifest. 

Specifically we consider the Fermi-Hubbard chain with 
hamiltonian 

H= J2 a(x-y)a(:z:)*.a(y) + ^^(a(:rr.a(x))' (1) 

with a{x)* ■ a{x) = at(^)* '^t(^) + '^li^)* o^li^)- ctix) is 
the hopping amplitude, satisfying a{x) — a{x)* , a{x) — 
a{—x), and A is the strength of the on-site interaction. 

H is integrable for the nearest neighbor hopping am- 
plitude, i.e. a(±l) = 1, a{x) — otherwise. For 



longer range hoppings H is commonly expected to be 
non-integrable. On the kinetic level, changing a amounts 
to changing the dispersion relation. Otherwise the struc- 
ture of the transport equation is not altered. Thus the 
issue of non-integrability is fairly accessible to the Boltz- 
mann kinetic equation. 

One aspect was studied in detail already in [6], where 
it was noted that for nearest neighbor coupling the 
Hubbard-Boltzmann equation has a much larger set of 
stationary solutions than usually anticipated. On the 
other hand for the domain of attraction of a non-thermal 
stationary state, the usual kinetic picture is valid. En- 
tropy is strictly increasing and the steady state is ap- 
proached exponentially fast. (We always work in the 
spatially homogeneous set-up.) Our goal here is to study 
the approach to stationarity once a is no longer of near- 
est neighbor type. As in [6] we will rely on numerical 
solutions of the Boltzmann-Hubbard equation and study 
two prototypical non nearest neighbor hoppings. 

(i) An additional next-nearest neighbor hopping term, 
i.e., a(0) = 1, a(±l) = -i, a{±2) = -f , and aix) = 
otherwise, with a tunable parameter 77 G K. The Fourier 
transform of a is the dispersion relation 

uj,j{k) = 1 — cos(27rfc) — 77Cos(47rfc). (2) 

The nearest neighbor case corresponds to 77 = 0. 

(ii) An exponential decay of higher-order hopping terms, 
i.e., q;(0) = —1, a{x) — — ie"'"'^' for x ^ 0, with a 
tunable parameter C > 0. The Fourier transform of a is 
the dispersion relation 
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uj^{k) = - ^ e"'^^' cos(27rjfc). (3) 
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The limit ( ^ oo corresponds to the nearest neighbor 
case after shifting and rescaling e^(l+u;^(fc)), while C ~^ 
allows for large hoppings of size l/C- 

Fig. 1 visualizes uj{k) for both cases i) and ii), as well 
as the "reference" nearest neighbor hopping model (black 
dashed line). Later the next-nearest neighbor hopping 
model is investigated numerically for a small rji — 2^ 
(dark green line in Fig. 1) as well as r] — ^ (hght green 
line) . 



w(k) 




FIG. 1. (Color online) The dispersion relation aj(fc) for the 
next-nearest neighbor model in Eq. (2) with rji = 2^ and 
'72 = I (green solid curves coinciding with the dashed line 
and with 2 local maxima, respectively), and for the exponen- 
tial hopping model in Eq. (3) with C — f (upper blue solid 
curve). All curves are shifted such that i^(0) = 0. The dashed 
curve shows the dispersion relation for the (reference) nearest 
neighbor hopping model. 

Our goal is to study the dynamics of the Hubbard chain 
at small interaction in dependence on ry, respectively C. 
For this purpose, in Section II we first recall the struc- 
ture of the corresponding Boltzmann transport equation. 
The collision rules for quasiparticles are implicitly deter- 
mined by conservation of momentum and energy, which 
will be discussed in Section III. The numerical scheme is 
explained in Section IV, which is the technical backbone 
of our investigations. We this tool we study the approach 
to a stationary state, see Section V, and its dependence 
on the collision rules, in other words on the dispersion 
relation. 

In [7] mass diffusion in dependence on 77 was studied for 
a "toy" linear transport equation. The divergence of this 
transport coefficient as 77 — ?> is related to our findings 
for the full Boltzmann equation. 



II. THE BOLTZMANN-HUBBARD EQUATION 

We briefly recall the structure of the Boltzmann- 
Hubbard equation, see [6] for details. For the Fourier 
transformation we use the convention 



Then the first Brillouin zone is the interval T = [—5,5] 
with periodic boundary conditions. The dispersion rela- 
tion Lu(k) — a(k) and, up to a constant, in Fourier space 
H can be written as 




+ ^ f d^k6{k)a^{kiy a^ika)* ai{k^)a^{ki) 

with k = fci + fc2 — fca — fc4 mod 1 and d*k = 
dki dk2 dk^ dk^. 

To arrive at the kinetic equation, we assume that the 
initial state of the chain is quasifree, gauge invariant, and 
invariant under spatial translations. It is thus completely 
characterized by the two-point function 

{dAky dr[k')) = 5{k ~ k')W„r{k). (6) 

It will be convenient to think of W{k) as a 2 x 2 ma- 
trix for each A: G T. Then, in general, W{ki)W{k2) ^ 
W{k2)W{ki) and every argument of standard kinetic the- 
ory has to be reworked. By the Fermi property we have 
< W{k) < 1 as a matrix for each k. In particular, W 
can be written as 

W{k)^ J2 e^{k)\k,a){k,a\, (7) 

where \k,a) for a € {t,i} is a /c-dependent basis in spin 
space and Sa are the eigenvalues with < £cr < 1- 

At some later time t the state is still gauge and trans- 
lation invariant, hence necessarily 

(a^(fc, t)* ar{k', t)) = d{k ~ k')W„r{k, t). (8) 

In general W{t) is a complicated object, but for small 
coupling A the quasi-free property persists over a time 
scale of order A~^, a structure which allows one to obtain 
the kinetic equation by second order time-dependent per- 
turbation theory. More details can be found, e.g., in [8- 
10]. Here we only write down the resulting Boltzmann 
equation 

^W{k,t)^C,[W]{k,t)+C^[W]{k,t)=C[W]{k,t), (9) 
at 

which has the structure of an evolution equation and 
has to be supplemented with the initial data W{k, 0) — 
W{k). 

The first term is of Vlasov type, 

C,[W]{k,t) = -i[H,s{k,t),W{k,t)], (10) 

where the effective hamiltonian Hcs{k, i) is a 2 x 2 matrix 
which itself depends on W . More explicitly, 

i?cff,l = j ^ dk2dk3dki Sik) V (i) 

X {W3W4-W2W3-W3W2-tT[W4]W3+tT[W2]W3 + W2) . 

(11) 
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Here and later on we use the shorthand W = 1 — W, 
Wi = W{ki,t), H,s,i = H,g{kut), w = w(fci) + w(fc2) - 
^(fcs) — uj{ki). Since is 2 x 2 matrix- valued, tr[-] is 
the trace in spin space. Finally V denotes the principal 
part. Since the fcs, integration can be interchanged, 
Hctf = H*g, as it should be. 

There are many different ways to write the collision 
term Cd- We choose a version which separates the various 
contributions into gain and loss term. Then 



Cd[Wh 



TT / dk2dk^dk4^5{k) 5(uj) 



{A[WU3i + A[W] 



1234; J 



(12) 



where the index 1234 means that the matrix ^[M^] de- 
pends on fci, ^2, fca, and k/^. Explicitly 

^[W']l234 = -WiW2Wi + Wa tr[VK2W^3] 

- {WiW:i - WaW2 - W2W3 + Wa tr[W2] 

- Wa trlWs] + tT[W3W2]}Wi (13) 

with the first two summands the gain term and {...}W^i 
the loss term. The gain term is always positive definite, 
as implied by the inequality 



A tr[BC] + C tr[BA\ - ABC - CBA > 



(14) 



valid for arbitrary positive definite matrices A,B,C. 
Thus if an eigenvalue of W{k,t) happens to vanish, the 
gain term pushes it back to values > 0. A similar argu- 
ment can be made for W{k,t), implying the propagation 
of the Fermi property [10], to say: if at t = one has 
< W{k) < 1, then the solution to (9) also satisfies 
< W{k,t) < 1. 
In general, " spin" , 



dkW{k,t) 



and energy 



dkuj{k)tr[W{k,t)] 



(15) 



(16) 



are conserved. In the long time, W{k, t) will become di- 
agonal in the conserved spin basis. Each component has a 
Fermi-Dirac distribution with common temperature and 
destined chemical potentials, which then is precisely in 
accordance with the parameters from the conservation 
laws. 

For the nearest neighbor model, one has the additional 
conservation law 



d_ 
dt 



(tr[W^(fc, t)] - tT[W{\ - k, t)]) = 0. 
All stationary states are necessarily of the form 



(17) 



(18) 



with f{k) = — /(^ — k). Wnth is an equilibrium state if 
/(fc)=/3w(fc). 

The entropy of the state W is then defined by 



S[W] = - dA:i(tr[M/ilogTyi] -|-tr[M^ilogT^i]). (19) 

in accordance with an ideal Fermi gas. It is easily checked 
that the entropy production cr > 0, 



a[w] = -S[W] 



Jd 



dki tr[(logiyi -logll"i)C[VF]i]. 



(20) 



The H-theorem asserts that 

a[W] > for ah W with 0<W <1. (21) 

III. COLLISIONS 
A. Next-nearest neighbor model 

The starting point is to investigate the kinematically 
allowed collisions 5{k)S{Lu^). Using momentum conserva- 
tion fc = mod 1 and defining S12 = /ci -I- ^2 = ^3 + ^4, 
Afci2 — ^(fci — k2) and Ak^A — 5(^3 ~ ^a), one arrives at 
the factorization 

(22) 
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FIG. 2. (Color online) Contour (green straight lines) and 
gradient (gray vectors) of the next-nearest neighbor energy 
conservation contour = (with rj = ^) for fixed ki — ^ 
and after eliminating k2 . The vertical and horizontal lines, 71 
and 72, are the contours ^3 = ki and fc4 = fei, respectively. 
The contour 7eiiip disappears when \r]\ < |. 
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with the factors 

^bas = 4sin(7r(fci - k^)) sin(7r(fci - k^)) 
= 2 (cos(27rAfc34) - cos(27rAfci2)) 



and 



:^add.r, = cos(7rsi2) +?7 COs(27rSi2) 

X (cos(27r Afci2) + cos(27r Afc34)) 



(23) 



(24) 



Eq. (22) is of similar form as [6, Eq. (34)], except for the 
additional 77-dependent term in Wg^^d rj- particular, the 
"trivial" solution paths fca = ki (denoted 71 ) and k^ = ki 
(denoted 72) remain unaffected by 77. A sign change of 
77, i.e., 77 — > —77, corresponds to fc^ — > fc^ + 4 since this 
transformation sends S12 -> S12 + 1 and cos(7rsi2) — >■ 
— cos(7r S12), while the other cosine terms in Eq. (24) are 
unaffected. Thus without loss of generality one may as- 
sume that 77 > 0. 
We decompose 

A[W] 

1234 1234 — -^quad[W^] 1234 +-4tr[W^] 1234 (25) 

with 



Auad[M^]l234 = -WiW^W^Wi - 
+ WiW3W2Wi 



WiW2W3Wi 

'W4W2W3W1, (26) 



Ar[W]l234 = {W1W3 + M^3W^l)tr[W^2M^4] 

- (W'iM^3 + W^3W^i)tr[T4^2W4]. (27) 

The discussion of the integration along 71 , 72 follows the 
same line as in [6]: .Aquad is zero along both 71, 72, but 
^tr is zero along 71 only. 

Concerning the factor uj^^^^ in Eq. (22), the con- 
tour w^dd.r) = splits into two parts, denoted 7diag 
and 7oiiip, see Fig. 2. Using the identity cos(27rsi2) — 
2cos^(7rsi2) 1 and solving for S12, one arrives at 



, , 1 / ±Vl + 2r2- 1 

Si2(f') = — arccos 

TT \ 2r 



for 7diag and 7oiiip, respectively, where 

r ^ Arj (cos(27r AA;i2) -I- cos(27r Afc34)) . 



(28) 



(29) 



The argument of the arccos function in Eq. (28) should 
be in the interval [—1,1], which is always satisfied for 
7diag- However, on 7cnip this constraint leads to the con- 
dition \r\ > 2. Thus we conclude that the contour 7ciiip 
disappears for I77I < 4 since by Eq. (29), |r| < 4 I77I 2 < 2. 

As a remark, Taylor-expansion at r = of Eq. (28) on 
7diag gives 



, , 1 r llH 



(30) 



In particular, we reobtain the constant S12 = ^ for the 
next-neighbor case 77 = 0. 



In summary, the contour 7diag is deformed as compared 
to the nearest neighbor case (compare with [6, Fig. 2]). 
The additional collision channel 7eiiip appears when \ri\ > 
4 . The gradient vector field of oj^ (gray vectors in Fig. 2) 
is noticeable different compared to the nearest neighbor 
case. 



B. Exponential hopping 



We analyze the kinematically allowed collisions 
6{k)d{uj_^) for the dispersion relation in Eq. (3). A short 
calculation shows that Eq. (3) can be written as 



UJ({k) 



sinh(C) 



cosh(^) — cos(27rfc) 



(31) 



Again using the momentum conservation fc = mod 1 
and some trigonometric identities, one arrives at the fac- 
torization 



= ^ sinh(C) CJbas ^add.C 

X j Y[ (cosh(^) — cos(27rfci)) 



(32) 



0.5 



- - - ' y 

. I - n 
y - ■ 


■ , , ^ . 
- V / >\ ' - 

■ rdiagV 

' / ■ x\ - 


\ ) ' — 
- • - - ' >^ 


. . . ^ / 

- \ » * * * t \ \ 1 

X ' • ■ - - 



0.5 

^3 



FIG. 3. (Color online) Contour (blue straight lines) and gra- 
dient (gray vectors) of the exponential decay energy conser- 
vation Lu^ — (with C = §) for fixed = |f and after 
eliminating ^2. The vertical and horizontal lines, 71 and 72, 
are the contours fcs = ki and ki = ki, respectively. Compare 
with Fig. 2 corresponding to the next-nearest neighbor case. 
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with the same factor Wj^j^g as in Eq. (23), and 

^add.C = ~ cos(7rsi2)^ + cos(7rsi2) 

X (l + cosh(C)^ + cos(27rAfci2) cos(27rAfc34)) 
— cosh(C) (cos(27rAfci2) + cos(27rAfc34)) . 

(33) 

^add ^ = is a cubic equation for cos(7rsi2), which can be 
solved analytically in closed form or numerically by a few 
Newton iteration steps. There is only a single real- valued 
solution, which we (again) denote by 7diag (the context 
will resolve any ambiguity to the next-neighbor case). 
Fig. 3 visualizes the contours iO(- = 0, which resemble the 
next-nearest neighbor case except that 7oiiip is missing 
and 7diag is slightly distorted. One notices that 7diag 
and 7i seem to intersect at (fc3,fc4) = (fci,0). This is 
no coincidence, since in the limit ^ — > 0, the equation 
^add c ~ admits a solution Afci2 = Afc34 = which 
is equivalent to fci = and fc2 = fc4 = 0. Similarly, 7diag 
and 72 intersect at fci = fc4, /c2 = ^3 = when — >■ 0. 
The nearest neighbor case [6] corresponds to the limit 
— >■ cx): namely, dividing Eq. (33) by cosh(C)^ (which 
leaves the solutions of lo^^^^ c ~ ^ invariant) and letting 
( — > cxD, only the term cos(7rsi2) remains. 



D. Integrable models 

The Hubbard chain is integrable for pure m-th neigh- 
bor hopping models with dispersion relation 



^m{k) — — cos(27rmfc). 



(36) 



Similar to nearest neighbor hopping (m = 1), the energy 
conservation factorizes as 



4 cos(7rm(fc3 -I- k^)) 

X sin(7rm(fci — ^3)) sin(7rm(A:i — k/C})- 



(37) 



Accordingly, the collision contours are re-scaled by the 
factor m. 

There is an infinite number of energy-like conservation 
laws: Let 5 : T -> M with g{k + ^) = g{k) for all A: e T, 
as well as g{k) = —g{^ — A). Then 



dt 



dkg{k)tT[Wik)] = 0, 



(38) 



which follows by an appropriate interchange of the inte- 
gration variables fci, . . . , ^4. Note that g{k) is completely 
determined by prescribing g{k) for fc e [ 



1 1 

4m ' 4m J 



IV. NUMERICAL PROCEDURE 



C. Stationary solutions 



A. Contour integrals of the dissipative collision 
operator 



In the spatially homogeneous case, the conventional 
wisdom is that the stationary solutions of the kinetic 
equation coincide with thermal equilibrium. This should 
hold also if in (1) the lattice Z is replaced by the d- 
dimensional lattice Z''. As proved in [6], for a general 
dispersion relation and in arbitrary dimension the prob- 
lem of classifying all stationary solutions can be reduced 
to finding the set of all collision invariants, i.e., solutions 
to 



$(fcl) $(fc2) = *(fc3) + *(fc4) 



(34) 



on the manifold {{ki, ^2, k^, ^4) | fc = mod 1, ij — 0}. 
The obvious solution reads 



^k) = ^{uik)-^l^), 



(35) 



which corresponds to thermal equilibrium. Thus the is- 
sue reduces to whether there are further collision invari- 
ants. For dimension d > 2 a proof under fairly general 
conditions is available [11]. For d ^ 1, there could be too 
few collision channels to reach thermal equilibrium. An 
example is the nearest neighbor Hubbard chain. There 
is then no 7onip and 7diag is linear. As a consequence 
additional collision invariants can be found. Our numer- 
ical simulations indicate that a slight curvature of 7diag 
suffices to limit the set of collision invariants to the ones 
listed in (35). 



The following discussion applies to both the next- 
nearest neighbor and exponential model. Ideally, the nu- 
merical discretization of the integration contours (Fig. 2 
and 3) should preserve the spin and energy conserva- 
tion laws. These conservation laws result from the inter- 
changeability /ci o ^2, ^3 O k^ and the pairs {fci, A:2} -H- 
{fc3, fc4}. For the contours 71 and 72, we can proceed as 
in [6] using a uniform grid for the k variables. However, 
the contours 7diag and 7oUip require more sophistication: 
to adopt the symmetries in the numerical discretization, 
we first rewrite the dissipative collision evaluated at k: 



TT / dkidk2dk-idki 5{k) (5(w) 5{ki - k) {A[W] -f A[W]*) 

= TT / dAki2dAk34 / (isi2 ds34 <5(si2 - S34) '5(w) 

X ,5 (S12/2 + A/C12 - k) {A[W] + A[W]*) 

= TT dAki2 (iAfc34 / dsi2 S{ui) 
Jt^ Jt 

X S (S12/2 + Aki2 - k) {A[W] + A[W]*) , 



(39) 



where we have used the substitution 

1 



S12 = ki + k2, Afci2 ^ fc2), (40) 

534 = ^3 + ^4, Afc34 = ^(fc3 - fc4). (41) 
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In the following we are only concerned with the integra- 
tion along the contour 7ciiag or 7eUip- The S12 integral 
can be eliminated the via namely, 



dsi2(5(w) = |asi2Wp\ (42) 



Thus the last integral in Eq. (39) becomes 

TT / (iAfci2dAfc34 |(9si2a;|~^ 
Jt2 (43) 

X 6 (S12/2 + Afci2 - k) {A[W] + A[W]*) , 

where S12 depends on Afci2 and Ak^^ via Eq. (28) or 
^add.c = in Eq. (33), respectively. Numerically, we 
discretize the integral in (43) by a uniform grid: 



Afci2 = 



J 



n n 
2'^2 



1, 



- 1 



(44) 



(same for Ak^^) with fixed n = 128 in our case. Note that 
fci O k2 corresponds to Afci2 O — AA:i2 and likewise 
for Afc34, and that {ki,k2} -H- {k^^k^} corresponds to 

Afci2 O Afc34. 

So far we have not taken the (5-function in Eq. (43) 
into account, for which we use the following approach: 
we want to determine the cumulative contribution to the 
collision operator at the uniform fc-grid points k — ^, 
j = 0, l,...,n — 1. We do not resolve the (5-function 
exactly; instead, for each term 



A^TT \ds12 wf^ {A[W] + A[WY 



(45) 



evaluated at discretized Afci2, Afc34, we first choose k 
i such that 



k < S12/2 + Afci2 < k 



1 



(46) 



Then we add v^A to Cd[W]{k) and (1 - v)^^A to 
Cd[VF](fc + -), with G M chosen such that 



a;(si2/2 + Afci2) = ly uj{k) + {l-v)uj [k + - j . (47) 

By this approach, the numerical scheme preserves the 
spin and energy conservation laws. 

In summary, our numerical method approximates 
Cd[W](fc) (and thus W{k,t) for the next time step) at 
the uniform fc-grid points k — ^. However, the discretiza- 
tion (44) of the terms in Eq. (45) requires evaluation of 
W{k) at isi2 ± Afci2 and ^342 ± Afc34, which are (in 
general) no uniform grid points ^. We solve this issue 
by polynomial interpolation of order 3 (precomputing di- 
vided differences based on W{k) at fc = ^). 



B. Mollifying the collision operators 



We use the same mollification scheme as in [6] to avoid 
the infinities resulting from |9si2 a;| in Eq. (45) and the 



principal value of 1/w in the effective hamiltonian (11) 
Concretely, we replace 



-1/2 



\dsi20o\-^^ {\dsi2Ui\^ + e') ■ , (48) 
and for the conservative collision operator 



V 



(49) 



with finite e > 0. In our case, we use e = for the 
simulations in section V. Note that Eq. (49) becomes an 
exact identity when taking the limit e — > 0. 

While the mollification parameter e is required to avoid 
infinities, it has to be chosen somewhat arbitrarily. We 
briefly quantify the effect of different values ei = 
^2 — jq and £3 = ^ in Fig. 4. The curves show the 
Hilbert-Schmidt difference \\W{t) - W{0)\\ between the 
current and an initial Wigner state (see section V) up to 
t — 2 (next-nearest neighbor model with rj = The ef- 
fects of different mollifications are quite noticeable for the 
time interval shown. On the other hand, the curves ap- 
proach each other for larger t since all Wigner states even- 
tually converge to the same thermal equilibrium state. 
Thus it is reasonable that the asymptotic convergence to 
equilibrium hardly depends on the mollification. 

l|W(t)-W(0)|| 




FIG. 4. Effect of different mollification parameters ei = 
(upper dark gray curve, used for the simulations in section V) , 
e2 = jfj (middle curve) and 63 = | (lower light gray curve). 
The difference to the initial W{k, 0) is quantified by the 
Hilbert-Schmidt norm. 



C. Solving the Boltzmann equation 

Departing from [6] , we avoid the Strang splitting tech- 
nique for treating Cd and Cc separately, but simply use 
the explicit midpoint rule for C = Cd+Cc- As advantage, 
this approach exactly preserves the spin and energy con- 
servation laws. The more laborious time evolution step 
for Cc in [6] did not show any noticeable differences. 
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D. Implementation details 

We have implemented the numerical scheme described 
so far in plain C code, with a custom struct for complex 
Hermitian 2x2 matrices (with 4 double values for the 
real diagonal entries and the complex 1,2 entry). The 
implementation is designed such that the intermediate 
steps always deal with Hermitian matrices. For example, 
the commutator S] and anticommutator {A, B} for 
Hermitian A, B is again Hermitian and can directly be 
calculated from the matrix entries of A and B, without 
resorting to the products AB or BA. Similarly, a custom 
function calculates the sum of triple products ABC + 
CBA directly from the matrix entries, which is again 
Hermitian when A, B,C are. 

We use the MathLink interface to make the numerical 
procedures conveniently accessible from Mathematica. 

The C implementation comes with a noticeable per- 
formance increase: on the same hardware as in [6] (Intel 
Core i7-740QM Processor, 6M cache, 1.73 GHz), a sim- 
ulation run with the same parameters as in [6] now only 
takes several seconds, as compared to 6 h for the Mathe- 
matica implementation in [6]. 



V. SIMULATION 
A. Initial Wigner state 

Our goal is to investigate the effects of the different dis- 
persion relations uj{k) in Fig. 1. We start with a (rather 
arbitrary) initial condition W{k, 0) shown in Fig. 5. The 
bright and dark cyan lines represents the real diagonals, 
and the dark and light red oscillatory functions the real 
and imaginary part of the off-diagonal |t)(4-| entry, re- 
spectively. The eigenvalues of W{k, 0) are in the interval 
[0, 1] for each A; G T, as required by the Fermi prop- 
erty. W{k, 0) is continuous on T. The analytic formula 
of W{k,0) can be found in appendix A. 



W(k,0) 

1 r 




FIG. 5. (Color online) The initial state W{k, 0) used for all 
simulations in this section. The cyan (upper) curves show the 
real diagonal entries, and the darker and lighter red curves 
the real and imaginary parts of the off-diagonal |t)(i| entry, 
respectively. 



As illustration. Fig. 6 visualizes the 3-dimensional 
shape of the collision manifolds 7diag and 7oiiip for the 
next-nearest neighbor model with rj — ^. (Note that 
Fig. 2 is the intersection of Fig. 6 with the hyperplane 
ki = ||.) To illuminate the effect of the dissipative col- 
lision operator Ca, the colors in Fig. 6 encode the Bloch 
vector of A[W] + A[W]* for the initial state W{k,0), 
where the red, green and blue colors correspond to the 
X, y and z components of the Bloch vector, respectively. 








FIG. 6. (Color online) 3D shape of the 7diag and 7ciiip collision 
manifolds for the next-nearest neighbor model with f? = | (as 
in Fig. 2). Color encodes the Bloch vector of A[W] + A[W]* 
for the state W{k, 0) shown in Fig. 5. 



B. Stationary states 

For the given initial W{k, 0), one can obtain the corre- 
sponding stationary state (which is different for different 
dispersion relations) from the conservation laws Eq. (15), 
(16) and (17). We will discuss four different models ac- 
cording to Fig. I: the nearest neighbor case (77 = 0), the 
next-nearest neighbor model with a small perturbation 
771 = 2^ and a larger ?72 = ^ (such that the 7oiiip colli- 
sion path opens up), as well as the exponential hopping 
model with C = §■ The corresponding stationary states 
are distinct. 

For the nearest neighbor model, the stationary solution 
is a non-thermal state of the form 

VF„th(fc)= Yl (c^''^-'^' + l)"'k)(fT|, (50) 
i.e., a real diagonal /c-dependent matrix, where the func- 
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(c) thermal equilibrium state for next-nearest 
neighbor hopping 
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(d) thermal equilibrium state for next-nearest 
neighbor hopping 



FIG. 7. (Color online) Diagonal matrix entries of the stationary states corresponding to the initial W{k,0) in Fig. 5, for the 
nearest neighbor hopping model (a), for the exponential hopping model with C = | (b), and for the next-nearest neighbor 
model with rj — ^^jg (c) and ?7 = | (d). The off-diagonal matrix entries are zero. 



tion / satisfies the symmetry property f{k) = — /(^ — k)- 
f is obtained mimerically, and the corresponding Wnth(^) 
shown in Fig. 7a. The cyan Unes visuahze the diagonal 
entries. 

With a perturbation rj =/= O'm the next-nearest neighbor 
model, the stationary states become thermal states of the 
form 

WtKnik)= E (e^^"''('=)-^"' + l)"V>('^l- (51) 

Figs. 7c and 7d visualize Wth,))(fc) calculated from the 
initial W{k, 0). Compared to f{k), the term /3 w^(fc) lacks 
the symmetry f{k) + f{^ — k) — 0. Note that even for 
77 0, in general W^th,jj(^) does not converge to Wnth- 

The numerically obtained values of /3 and /Xo- in 
Eq. (51) are summarized in the following table: 









= 0.005 


0.650 0.949 0.061 


m 


= 0.5 


0.752 0.972 0.176 



The stationary state of the exponential hopping model 



is a thermal equilibrium state of the form 

mh,c(fc)= E (e^^"^^')"^'^ + l)"V)(^l (52) 

as shown in Fig. 7b. The corresponding parameters are 
/3 — 1.00, fj,^ — —1.00 and /i^ = —1.60. The peak around 
k = becomes sharper when ^ decreases. 

C. Exponential convergence, fast and slow motion 

We pick the entropy as representative measure of con- 
vergence to stationarity. In our numerical simulations, 
we observe exponential convergence, i.e., the entropy dif- 
ference 

S[Wst] - S[W{t)] ~ e-'^* (53) 

for large times t, where Wst denotes the respective sta- 
tionary state. The following table summarizes the decay 
rates k obtained from a least-squares fit in logarithmic 
representation: 
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S[W(t)], ?7=0.005 
1.297 



1.2966 



1.2962 
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(a) entropy increase 
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|W(t)i2l, ?7=0.005 
0.10 



0.08 



0.06 
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0.02 



10 20 30 40 50 

(b) convergence of the ofE-diagonal entries 



FIG. 8. (Color online) Entropy increase for the next-nearest neighbor model with small 77 



(green curve). The red curve 



shows the entropy of the corresponding equilibrium state, and the dashed black curve the entropy of the stationary nearest 
neighbor state. The entropy grows very slowly after t ~ 10. 





nearest 


771 = 0.005 


772 = 0.5 


C = 0.4 




0.852 


0.001 


0.0676 


0.0530 



One notices that the convergence rate is highest for the 
nearest neighbor model, and lowest for the next-nearest 
neighbor model with small 771 = 250- We investigate the 
latter case in more detail. The green line in Fig. 8a shows 
a closeup of the entropy in dependence of t, and 

the red line the entropy value S[Wth,7j{k)] = 1.297 of 
the corresponding stationary state. For comparison, the 
dashed line is the entropy of the non-thermal stationary 
state (r/ = 0). One notices that the entropy grows much 
faster when t < 10 and then reaches a plateau, where 
it approaches the asymptotic red line very slowly. This 
observation suggests the following dynamical picture: In 
the phase space for (9) there is the slow manifold con- 
sisting of Wigner functions of the form (50). A general 
initial state, W, will rapidly move towards the slow man- 
ifold, and will arrive there at a state W{t^), where in 
general W{t^) ^ Wnth- From there on there is an effec- 
tive dynamics on the slow manifold with initial Wigner 
function W{t^.). This can be seen in Fig. 8b. To obtain 
the evolution equation in the slow manifold, we treat the 
off-diagonal entries as small perturbation, 

Wik) = W^{k) + SW'^^ik) (54) 

with < (5 < 1, Ty°(fc) the diagonal part and W^°{k) 
the off-diagonal part. The effective dynamics for the state 
is driven by 



C[W^ +SW°^]{k). 



(55) 



Since the conservative collision operator Cc in Eq. (10) is 
defined by a commutator, it holds that 



C,[W]{k,t) = 0{S). 



(56) 



Thus the Boltzmann differential equation is to zero-th 
order in S governed by the dissipative part coupling the 
tt and IX correlation functions: 



where 



C^[W^^](fc,t) = TT / dk2dk3dk4S{k)S{uj_) 
Jt3 



dt 



W^'ik^t) = C^[W'']{k,t) + 0{S), 



(57) 



(58) 

The differential equation for W^^ is given by interchang- 
ing tt and II in Eq. (58). We suspect that this is the 
effective equation for the slow-motion dynamics. 

The concept of different dynamical regimes is sup- 
ported by Fig. 9: The dark gray points represent the 
inverse asymptotic decay rates 1/k for the next-nearest 
neighbor model in dependence of 77. For comparison, the 
light gray points show the initial decay rates at W{k, 0). 
One observes that initial and asymptotic decay rates are 
clearly separated. 



VI. CONCLUSIONS 

On the level of the Boltzmann-Hubbard equation one 
can easily destroy integrability by going beyond the next- 
nearest neighbor hopping. The structure of the kinetic 
equation is not touched, but through modifying w one 
changes the set of allowed collisions. The consequences 
on the dynamics are in accordance with text book wis- 
dom. In the integrable case the collision rule has a high 
symmetry and, while there is still exponential conver- 
gence and non-zero entropy production, in general one 
reaches a nonthermal state of the form (50). Any tiny 
modification of lo restores the physically expected ther- 
malization to the Fermi-Dirac diagonal Wigner function. 
For large modifications we find again exponential fast 
convergence. However, for a small perturbation of w, 
we clearly demonstrated two time scales, a rapid con- 
vergence to quasi-stationarity and a slow convergence to 
thermal equilibrium. 

Our model is fairly simple, but serves as an example 
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(a) exponential decay rate in dependence of r) 
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FIG. 9. (Color online) (a) Inverse exponential decay rate 1/k of the entropy diflerence in dependence of r; (next-nearest neighbor 
model), obtained from a least squares fit as exemplified in (b). The upper dark gray points in (a) correspond to the asymptotic 
decay rate for large t (dark dot-dashed line in (b)), and the lower light gray points to the initial decay rate at i = (light 
dashed line in (b)). 



where the approach (and non-approach) to thermal equi- 
librium can be studied in detail. 



Appendix A: Analytic formula of W{k, 0) 

For the sake of reproduceability, the analytic formula 
of the initial Wigner state W{k^ 0) used in the simulations 
(section VA) reads as follows: 

+ ^ ^18 cos (7r(6fc + 1/7)) - 14cos (67r(fc - 1/7)) 
+ 27 (cosh(l) - cos(47rfc))~^ (e"^/^ cos(27rA:) 

+ cos(47rfc) - e^/^ cos(67rfc) + e"^^ ^ , (Al) 



W^^{k, 0) = ^ (9 sin(e8''^'=) - (1 + i) cos (67r(fc 1/7)) 

+ sin (7r(3fc + 1/14)) sin (37r(fc - 1/7))) (A2) 
together with 

W^^{k,Q)^Wn{k,0)\ (A3) 

and 

w^u(fc,o) 

^ ^gi(cos(27rfc)-cosh(2/5))-^ sinh(2/5)+ii + ij'^ 

-I- 1^ ^14 COS (67r(fc - 1/7)) - 18 COS (7r(6 A: + 1/7)) 
-I- 27 (cosh(3/2) - cos(47rfc))"^ (^e""/i° cos(27r/c) 

+ cos(47rA:) - e^/^ cos(67r/c) - e"^/^) ^ . (A4) 
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